% This function takes estimates and imputed networks and imputes outcomes
 
function [y1_imp] = Fn_ImputePE_Mar2022(PE_out, A_hat, W_ij, Cons, y1_lag)

% Setup
    LMH_educ = [Cons.L_educ, Cons.M_educ, Cons.H_educ];
    X_0 = [ones(Cons.n,1), Cons.P, W_ij*Cons.P];
    X_0_educ_LMH = [];
    for z=1:3
        X_0_educ_LMH = [X_0_educ_LMH, bsxfun(@times,LMH_educ, X_0(:,z))];
    end

    LMH_gender = [Cons.L_gender, Cons.M_gender, Cons.H_gender];
    X_0 = [ones(Cons.n,1), Cons.P, W_ij*Cons.P];
    X_0_gender_LMH = [];
    for z=1:3
        X_0_gender_LMH = [X_0_gender_LMH, bsxfun(@times,LMH_gender, X_0(:,z))];
    end
    
    
%%%%% Without A's
    X_educ_m1 = X_0_educ_LMH;
    y1_educ_imp_m1 = X_educ_m1*PE_out.alpha_educ_m1 + sqrt(PE_out.Sigma2_educ_m1)*randn(Cons.n,1);
    y1_imp.y1_educ_m1 = Cons.y1_educ_MISS.*y1_educ_imp_m1 + (ones(Cons.n,1)-Cons.y1_educ_MISS).*y1_lag.y1_educ_m1;

    X_educ_m3 = [X_0_educ_LMH, bsxfun(@times,LMH_educ,Cons.y0_educ), bsxfun(@times,LMH_educ,W_ij*Cons.y0_educ)];
    y1_educ_imp_m3 = X_educ_m3*PE_out.alpha_educ_m3+sqrt(PE_out.Sigma2_educ_m3)*randn(Cons.n,1);
    y1_imp.y1_educ_m3 = Cons.y1_educ_MISS.*y1_educ_imp_m3 + (ones(Cons.n,1)-Cons.y1_educ_MISS).*y1_lag.y1_educ_m3;
    

    X_gender_m1 = X_0_gender_LMH;
    y1_gender_imp_m1 = X_gender_m1*PE_out.alpha_gender_m1 + sqrt(PE_out.Sigma2_gender_m1)*randn(Cons.n,1);
    y1_imp.y1_gender_m1 = Cons.y1_gender_MISS.*y1_gender_imp_m1 + (ones(Cons.n,1)-Cons.y1_gender_MISS).*y1_lag.y1_gender_m1;
 
    X_gender_m3 = [X_0_gender_LMH, bsxfun(@times,LMH_gender,Cons.y0_gender), bsxfun(@times,LMH_gender,W_ij*Cons.y0_gender)];
    y1_gender_imp_m3 = X_gender_m3*PE_out.alpha_gender_m3+sqrt(PE_out.Sigma2_gender_m3)*randn(Cons.n,1);
    y1_imp.y1_gender_m3 = Cons.y1_gender_MISS.*y1_gender_imp_m3 + (ones(Cons.n,1)-Cons.y1_gender_MISS).*y1_lag.y1_gender_m3;
    

    
    
%%%%% With A's
if var(A_hat)>1e-4    
    X_educ_m2 = [X_0_educ_LMH, bsxfun(@times,LMH_educ,A_hat), bsxfun(@times,LMH_educ,W_ij*A_hat)];
    y1_educ_imp_m2 = X_educ_m2*PE_out.alpha_educ_m2+sqrt(PE_out.Sigma2_educ_m2)*randn(Cons.n,1);
    y1_imp.y1_educ_m2 = Cons.y1_educ_MISS.*y1_educ_imp_m2 + (ones(Cons.n,1)-Cons.y1_educ_MISS).*y1_lag.y1_educ_m2;

    X_educ_m4 = [X_0_educ_LMH, bsxfun(@times,LMH_educ,A_hat), bsxfun(@times,LMH_educ,W_ij*A_hat), bsxfun(@times,LMH_educ,Cons.y0_educ), bsxfun(@times,LMH_educ,W_ij*Cons.y0_educ)];
    y1_educ_imp_m4 = X_educ_m4*PE_out.alpha_educ_m4+sqrt(PE_out.Sigma2_educ_m4)*randn(Cons.n,1);
    y1_imp.y1_educ_m4 = Cons.y1_educ_MISS.*y1_educ_imp_m4 + (ones(Cons.n,1)-Cons.y1_educ_MISS).*y1_lag.y1_educ_m4;
    
    X_gender_m2 = [X_0_gender_LMH, bsxfun(@times,LMH_gender,A_hat), bsxfun(@times,LMH_gender,W_ij*A_hat)];
    y1_gender_imp_m2 = X_gender_m2*PE_out.alpha_gender_m2+sqrt(PE_out.Sigma2_gender_m2)*randn(Cons.n,1);
    y1_imp.y1_gender_m2 = Cons.y1_gender_MISS.*y1_gender_imp_m2 + (ones(Cons.n,1)-Cons.y1_gender_MISS).*y1_lag.y1_gender_m2;
 
    X_gender_m4 = [X_0_gender_LMH, bsxfun(@times,LMH_gender,A_hat), bsxfun(@times,LMH_gender,W_ij*A_hat), bsxfun(@times,LMH_gender,Cons.y0_gender), bsxfun(@times,LMH_gender,W_ij*Cons.y0_gender)];
    y1_gender_imp_m4 = X_gender_m4*PE_out.alpha_gender_m4+sqrt(PE_out.Sigma2_gender_m4)*randn(Cons.n,1);
    y1_imp.y1_gender_m4 = Cons.y1_gender_MISS.*y1_gender_imp_m4 + (ones(Cons.n,1)-Cons.y1_gender_MISS).*y1_lag.y1_gender_m4;
    
    
end




 
end

